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The recently proposed Hysteretic Optimization (HO) procedure is applied to the ID Ising spin 
chain with long range interactions. To study its effectiveness , the quality of ground state energies 
found as a function of the distance dependence exponent, a, is assessed. It is found that the transition 
from an infinite-range to a long-range interaction at a = 0.5 is accompanied by a sharp decrease in the 
performance . The transition is signaled by a change in the scaling behavior of the average avalanche 
size observed during the hysteresis process. This indicates that HO requires the system to be infinite- 
range, with a high degree of interconnect ivity between variables leading to large avalanches, in 
order to function properly. An analysis of the way auto-correlations evolve during the optimization 
procedure confirm that the search of phase space is less efficient, with the system becoming effectively 
stuck in suboptimal configurations much earlier. These observations explain the poor performance 
that HO obtained for the Edwards- Anderson spin glass on finite- dimensional lattices, and suggest 
that its usefulness might be limited in many combinatorial optimization problems. 



I. INTRODUCTION 

The steadily increasing availability of powerful com- 
putational resources has allowed scientists and engineers 
alike to study ever more realistic and complex problems. 
Historically, a significant fraction of all available com- 
puter time has been used to traverse phase space in search 
of the optimal solution to a given problem. 

Particularly, in several physics domains [1], such as spin 
glasses, disordered materials and protein folding, one is 
interested in enumerating a large number of local minima 
of the energy landscape [2, 3], as that provides us with 
valuable information about its physical properties. Many 
algorithms and heuristics have been developed over the 
years to tackle this kind of problem. In this article we 
study a recently proposed algorithm known as Hysteretic 
Optimization (HO) [4-7]. HO is motivated by the physics 
of demagnetizing magnetic materials with a slowly os- 
cillating external field of decreasing amplitude. Similar 
to the thermal fluctuations in simulated annealing [8] or 
the activated dynamics of extremal optimization (EO) 
[9] , the drag created by the external field carries HO over 
energetic barriers. 

HO has proved efficient [6] at finding the ground state 
configurations of the well know Sherrington-Kirkpatrick 
(SK)[10] mean-field spin glass and of the classical Trav- 
eling Salesman Problem, but has, hitherto, been unsuc- 
cessful in searching the corresponding Edwards- Anderson 
spin glass on finite-dimensional lattices [11], as described 
in Chap. 10 of Ref. [1]. Similarly, it was shown in 
Ref. [12] that HO performs poorly for the random field 
Ising model (RFIM) on a one or three dimensional lat- 
tice. (The RFIM is a classic model for disordered mag- 
netic materials, but unlike the spin glass case there are 
polynomial-time algorithms to find global optima in the 
energy landscape, see Chap. 5 in Ref. [1].) Since efficient 
heuristics for hard (i. e. beyond polynomial) optimiza- 
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tion problems are still few and far between, especially 
for spin glasses, but also for many other combinatorial 
problems [13], promising new algorithms warrant careful 
investigation. Here, we explore the behavior of HO un- 
der variation of the search space characteristics that is 
representative of many problems. In particular, we ap- 
ply HO to a one-parameter family of spin glass problems 
[14-18] that interpolates between the characteristics of 
the SK model on one extreme and the EA model on the 
other. We observe that the break-down of HO for the 
EA model is intimately linked to the physics of intermit- 
tent events (i. e. avalanches) kicked-off by the external 
field in the hysteresis loop. We find a distinct cross-over 
between broadly-distributed avalanching dynamics in the 
SK-regime[19], connected with a high degree of intercon- 
nect ivity between variables and divergent energy scales, 
and sharply cut-off dynamics in the EA regime. Unfor- 
tunately, the need for strong interconnectivity between 
variables severely limits the applicability of HO with re- 
spect to combinatorial problems related to spin glasses of 
low degree such as satisfiability, partitioning, or coloring 
at their respective phase transitions [13]. 

This article is structured as follows, in Section II, we 
introduce the general Hysteretic Optimization procedure 
and apply it to a generalization spin glass model. In Sec- 
tion HI, we investigate in detail the avalanche dynamics 
during the hysteresis process to identify the reasons that 
lead to the breakdown of HO's performance. 



II. HYSTERETIC OPTIMIZATION 

For a magnetic material, such as an Ising system, to 
obtain a zero magnetization value, an a.c. demagnetiza- 
tion is performed. The sample is placed in an oscillating 
and slowly decaying magnetic field. As the amplitude of 
the external field approaches zero, so does the magnetiza- 
tion. At low enough temperature and slow driving, a dis- 
ordered systems gets dragged through a sequence of local 
energy minima. Based on this observation, Zarand et al 
[4] proposed Hysteretic Optimization (HO) as a general- 
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Figure 1: Simple sketch of a one- dimensional slice through 
the energy landscape for a generic optimization problem as in 
Eq. (1) without external field (bottom) and in Eq. (2) with 
external field (top). 




Figure 2: Systems described by the long-range spin glass 
model as a function of a. For a = we re-obtain the SK 
model and as cr ^ oo all the long range links become essen- 
tially negligible leaving us only with nearest-neighbor inter- 
actions. 



purpose local search heuristic [20] to explore the phase 
space of many combinatorial optimization problems. As 
an example, their study implemented the HO algorithm 
as listed in Tab. 1 

Our study here is focused only on finding the ground 
state (T = 0) energies of spin glasses, for which case we 
describe the implementation of HO in detail. 

In an Ising spin glass, each spin Si G {±1} is assumed 
to have a random bond Jij with other spins Sj: 

H = -^JijSiSj, (1) 

where the summation is taken over all pairs of spins. To 
find ground states of this spin glass with HO, we couple 
each spin cr^ to an external field of amplitude H. with 
a random sign G {±1}, which may be adjusted even 
during a single demagnetization run. The Hamiltonian 
of this extended system is then: 

nHo = n + HY,^iSi (2) 

i 

Physically, the second term in Eq. (2) distorts the energy 
landscape as shown in Fig. 1, allowing the system to 
escape local minima. Fluctuations due to the coupling to 
the external field can compensate for unsatisfied bonds. 
By varying the external field H and the random couplings 
, HO can force the system to explore a vast area of phase 
space in search of the optimal solution. 



Algorithm 1 Hysteretic Optimization 

1. Set H — Hi large enough such that Si — \f i. Set 

^{min} = ^ (= Hho\h^q)- 

2. Decrease H until one spin becomes unstable and allow 
the system to relax. liTi < ^{min}, set ^{min} = 'H. 

3. Optional: When H passes zero, randomize ^i, leaving 
the current configuration stable. 

4. At each turning point H = Hn = — 7n-ii^n-i, for < 
7n < 1, reverse the direction of H. 

5. Terminate when amplitude \Hn\ < i^{min}r 

6. Restart at 1 for Nmn times with a new, random set of 

7. Return the best ^{min} over all runs. 



Following the prescription of Algorithm 1 for the vari- 
ation of the external field each run, in effect, starts 
by exploring a large region of phase space which subse- 
quently decreases slowly. By varying the field between 
positive and negative amplitudes Hn = —jn-iHn-i, 
the runs repeatedly quench the system, following an ap- 
proach similar to the well known simulated annealing or 
tempering algorithms [8, 21, 22]. 

HO operates at T = , thus there are no thermal fluc- 
tuations and we can simply calculate the field necessary 
to make the next spin unstable and increase it to that 
value (within the j'^Hq limit) . Typically, HO is run 
with multiple restarts from the largest amplitude to in- 
crease the chances of finding a better approximation to 
the global minimum. Note, however, that each run itself 
has no stochastic element once the couplings to the 
external field are fixed. It is therefore useful to restart 
the demagnetization process repeatedly with a fresh set 
of random field directions (see item 6. in Algorithm 1). 
In fact, it is also possible to refresh the each time the 
external field H passes through zero during each run (see 
item 3 in Alg. 1). 

This algorithm has been very successful in determining 
the ground state energies of the Sherrington-Kirkpatrick 
spin glass [6] and reasonably efficient for the Traveling 
Salesman Problem[23], but there are few attempts to ap- 
ply it to other problems[12]. In this article we focus on a 
Ising spin glass on a one-dimensional ring with power-law 
interactions [14-18] defined by the Hamiltonian in Eq. (2) 
with bonds of the form: 

= ^- (3) 

Here, Cij are random variables drawn independently from 
a Gaussian distribution of zero mean and unit variance 
and 

is the distance between each pair of spins on the ring. By 
varying a we can interpolate between the all-to-all SK 
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Figure 3: Phase diagram, after [14-18]. 
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Figure 4: Percentage difference between the ground state 
found by Hysteretic Optimization and the one found by Ex- 
tremal Optimization. 

limit (cr = 0) and the nearest-neighbor EA limit (large 
cr), as shown in Fig. 2. This model has been extremely 
useful in elucidating the connection between mean-field 
and finite-dimensional spin glasses[17, 18]. 

As it has been shown in the literature [14, 24, 25], as a 
is increased this spin glass goes through several distinct 
phases, see Fig. 3. For < a < 0.5 the system is ef- 
fectively Infinite Range (IR). For all a, the singular part 
of the mean field transition temperature, Tj^^, is of the 
order [17].: 

N N 
i=2 i=2 

where [•]^^ denotes an average over disorder with [ef^ ]at; = 
1. This temperature becomes finite in the thermody- 
namic limit at a = 0.5, signaling a transition to a Long 
Range (LR) regime, where each node is able to see only 
a finite fraction of the rest of the system. At a = 1.0, Tc 
becomes zero , but the LR character is preserved until 
cr = 1.5. From this point on, the structure of the system 
is purely Short Range (SR), and each spin is connected 
only to 0(1) neighbors 

We have performed a benchmark study of the per- 
formance of HO on this spin glass model over a range 
of cr— values. The point of this study is not so much 
to tweak HO for optimal performance, but to obtain a 
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Figure 5: Average avalanche size as a function of a for differ- 
ent system sizes. 



clear assessment of its behavior under variation of this 
parameter. To this end, we generated a benchmark of 
instances of system sizes N = 32, 64, 128, and 256, for 
which we have obtained extremely good approximation to 
the ground state energy by alternate means. In this case, 
we have used the Extremal Optimization heuristic {EO) 
[9, 26, 27] but expanding a large amount of CPU time to 
ensure accuracy. In fact, using the implementation de- 
scribed in Ref. [9], EO has proven itself equally capable 
of approximating ground states in the SK model [27] as for 
the EA[9], and it appears to be much less dependent on 
cr. Although a direct comparison is not justified here due 
to the disproportionate run times used for EO, we have 
found that even at much more extensive runs, HO was 
not able to find the exact-known ground state of a one- 
dimensional spin glass with more than ^ 10^ spins[12].) 
In this set of instances, we have applied HO with a min- 
imal set of control parameters. We set 7 = 0.99 and 
for each instance in our set, we performed 10 different 
runs, each with a separate sets of that were kept con- 
stant throughout the entire run. The ground state en- 
ergy was taken to be the best value seen over 10 different 
quenches. This value was then averaged over 1000 differ- 
ent instances and compared with the results obtained by 
EO for exactly the same set of instances. 

In Fig. 4 we plot the quality of the solution obtained by 
HO as a function of cr. The "Error" is defined as the per- 
centage difference in ground state energy of the solutions 
found using HO relative to EO. Generally, the quality 
of the results found by HO diminishes for increasing sys- 
tem size for all cr, as can be expected with the limited 
CPU time (linear in N) apportioned to these runs. More 
noticeable is the ever more pronounced rise in error for 
cr > 0.5. To understand the physical reasons behind this 
behavior, we proceed to studying the dynamics of the 
system in the next section. 
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III. AVALANCHES AND CORRELATIONS 

Unlike the comparison in the previous section, we now 
focus exclusively on the intrinsic behavior of HO itself. 
We will pinpoint the causes of i^Os breakdown using 
quantitative measures. Using essentially the same pro- 
gram as previously, for each value of A/", we perform 2 
(undamped) hysteresis cycles each, but for a much larger 
set of 4-10^-y32/]V instances. Throughout, we set ^ = 1. 

When the hysteresis procedure described causes a spin 
to become unstable and get flipped, this may cause sev- 
eral other spins to become unstable, thus initiating an 
avalanche in the system. Each avalanche can involve a 
significant fraction of the number of spins in the system, 
including, on occasion, several flips of the same spin in 
a form of long-range self interaction. Avalanches have a 
wide range of sizes and can, in principle, be larger than 
the system size by flipping the same spin multiple times. 

As a first step in our analysis, we measure (5'(cr)), 
the average avalanche size as a function of cr at different 
system sizes. This measurement will help us determine 
what range of a we need to study, since it should become 
system size independent in the nearest-neighbor limit. As 
we show on the right hand side of Fig. 5, this happens 
near a ^ 2.0, thus restricting our interval of interest to 
a G [0,2], as expected from the literature [17]. 

We find that this quantity obeys an empirical scaling 
relation of the form: 

log (TV) 

where A (a), B (a) are linear functions of a, and / (A^, a) 
is plotted in Fig. 6. The overlap of all the curves in 
the interval a G [0, 0.5] means that all the N dependence 
has been captured by the N/ log^ (A^) term. This scaling 
should be compared with the N/ \og{N) scaling found 
for this behavior for a = , i. e. the SK model, by Ref. 
[19]. In fact, the emergence of log-corrections makes any 
definite determination of scaling behavior impossible over 
the range of system sizes N accessible here, and any of 
the following scaling relations should be viewed as purely 
phenomenological. 

The scaling becomes increasingly worse with a > 0.5, 
signaling a new N dependence. 

We believe this change in behavior at a = 0.5 is due 
to the topological change, from IR to LR, that occurs 
at this point (see the discussion in the previous section). 
Even though each node in the LR regime is still con- 
nected to other spins at arbitrarily large distances, its 
possible influence is now limited to a fraction of the total 
number of variables in the system, resulting in smaller 
avalanches. The avalanche size effectively creates a limit 
on the length of the jumps in configuration space that 
the system is capable of performing, forcing a less than 
optimal sampling of phase space, and increasingly poorer 
results. 

Avalanche sizes are determined by the total number of 
spin flips that occur. If the same spin happens to flip 



several times, then it will be counted multiple times as 
well, but we can also count the number U of just which 
spins flip at least once. The ratio S/U of the avalanche 
size, iS, over the number of unique spins flipped, /7, gives 
us a measure of how important loops are in the dynamics 
of the system, a large ratio will indicate that perturba- 
tions spread throughout the system and keep returning to 
the same spin, while a number close to unity would mean 
that avalanches propagate in just one direction and never 
double back. 

On Fig. 7 we plot {S/U (cr) — 1) for different system 
sizes. We find that the ratio between the size of the 
avalanche and the number of unique spins flipped is al- 
ways small and becomes approximately system size inde- 
pendent in the Short Range phase. This confirms, once 
again, that in this region the sphere of influence of each 
spin is very small, being limited practically only to Near- 
est Neighbors. 

This quantity obeys a phenomenological scaling rela- 
tion of the form: 

^{N,a)-l^ = VN[g{N,a)-A\og{N)] (5) 

where A ^ 4 x 10~^is a small constant and g (A/", a) is 
shown in Fig. 8. The scaling collapse of the data is 
very good up to near a 1.5 where the system acquires 
a purely short range behavior. [Clearly, this collapse is 
purely phenomenological, as the log-correction in Eq. (5) 
would ultimately overwhelm the constant term.] 

Finally, we study how the algorithm approaches the 
final configuration, at ~ and m = 0, by looking at 
the auto-correlation function given by: 

{SfSl)-{Sf){Sl) (6) 

where the indices denote summation over all spins i. In 
Eq. (6), we measure the overlap between the final con- 
figuration and those obtained a number of r complete 
cycles backwards in the past at their H = 0-crossing. 
Intuitively, we expect that the configurations seen at the 
beginning of the procedure (large values of r) will be 



N=32 




0.5 1 1.5 2 

a 

Figure 6: Corrections to scaling of the average avalanche size, 
rescaled according to Eq. (4). 
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Figure 7: Ratio between the avalanche size and the number 
of spins that were flipped. 




Figure 8: Scaling function g{N,a), as defined in Eq. (5). 



completely unrelated to the final configuration (r = 0), 
resulting in a value near zero for this quantity. However, 
as the algorithm takes its course and approaches its con- 
clusion, so too must the configurations start approaching 
the final one, corresponding to a value close to 1. The 
way in which it varies from values near to values close 
to 1 gives us information about the way exploration of 
configuration space occurs. The longer the period dur- 
ing which the correlations are close to 0, the larger the 
volume explored, and the faster it gets close to zero, the 
earlier the system restricts itself to a given region, thus 
limiting the quality of the solution it is able to find. 

In Fig. 9 we plot this quantity for the case of = 256, 
averaged over 1000 different instances for each value of a 
and with 10 different runs per instance. For small values 
of cr, the plateau at low correlations is extended (lower 
solid black curve), followed by an increase towards the 
value of 1 near the final stages r ^ 0. As a increases, 
the auto-correlations increase within the plateau which 
itself shortens, and the tendency towards 1 becomes no- 
ticeable right from the onset (upper solid red curve for 
(7 = 2). This is a clear demonstration of the ideas ex- 
pressed earlier, that the volume of configuration space 
explored becomes smaller with the decrease in avalanche 
size corresponding to increasing a. 




T 

Figure 9: Correlation with the final solution as a function of 
the time in the past. 



IV. DISCUSSION 

The performance of the Hysteretic Optimization pro- 
cedure for spin glasses was analyzed for a spin glass 
model that interpolates between systems with highly con- 
nected variables in the mean field limit and sparsely 
connected variables in the nearest-neighbor lattice limit. 
HO is shown to be very fast, but the quality of it's so- 
lutions quickly start decaying for increasing values of 
the distance dependence exponent, a. An analysis of 
the avalanche dynamics occurring in these systems re- 
vealed that the failure of HO is due to the truncation of 
avalanche size, and hence a limited exploration of the en- 
ergy landscape, that occurs when the system is no longer 
in the Infinite Range phase. 

The analysis of the behavior of the auto-correlation 
function with a confirmed this idea by showing that HO 
becomes stuck in a limited region of configuration space 
increasingly earlier for larger values of the distance de- 
pendence exponent, a. 

HO, being dependent on avalanches for its local search, 
cannot continue to work when the avalanches are no 
longer large enough to facilitate large jumps in configura- 
tion space. Any attempt to use HO in a finite connected 
system, such as an Edwards- Anderson spin glass or many 
combinatorial optimization problems[13], is, therefore, 
inefficient (see Chap. 10 in Ref. [1] and Ref. [12]). Our 
attempts to simulate sparsely connected systems, in this 
case 3SAT[13] and EA spin glasses with ± J— bonds, with 
discrete bond weights proved particularly unsuccessful. 
In such system, all variables only possess a finite (and 
typically, small) range of local field states to take on. For 
instance, in such an EA spin glass in d = 3 dimensions, 
all spins have exactly 2d-\-l = 7 states. Thus, a hystere- 
sis loop has just 7 jumps between full up- and and full- 
down saturation. At each jump, a finite fraction of spins 
flip simultaneously due to degeneracies, but mostly in an 
uncorrelated manner dictated by their local environment. 
An open question remains that concerns problems defined 
on random graphs of finite connectivity but with a con- 
tinuous distribution of bond- weights, such as the Viana- 
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Bray spin glass [28] with a Gaussian bond distribution. 
Unlike for a lattice, the number of neighbors increases 
exponentially with distance so that every node is con- 
nected to every other node with ~ In(A^) steps, and even 
small correlations could quickly span the system. One 
might expect that there would be a crossover between 
the average connectivity and In (A/") separating broadly 
distributed avalanches from localized ones, which would 
be very weak. Hence, HO might still work reasonably 
well in those systems down to low connectivities for most 
practical system sizes. In fact, our preliminary studies of 
Gaussian spin glasses on 3-connected graphs showed only 
minor deterioration in HO compared to EO for increas- 
ing system sizes (up to TV = 1023). Yet, an independent 
comparison of HO to itself within a one-parameter family 
of models in the spirit of our approach here would require 
much more simulation for variable connectivity. 



These results also highlight one important ingredient 
for any efficient algorithm or heuristic: the ability to 
travel between very distant regions of configuration space 
without being impaired by the large energy barriers that 
make such jumps energetically or entropically unfavor- 
able. This ability is only within HOs reach for Infinitely 
Range systems. 



V. ACKNOWLEDGMENTS 

We would like to thank Helmut Katzgraber for inspir- 
ing discussions and the Division of Materials Research at 
the NSF for their support under grant #0312150 and the 
Emory University Research Council for seed funding. 



[1] A. Hartmann and H. Rieger, editors. New Optimization 
Algorithms in Physics. Springer, Berlin, 2004. 

[2] H. Frauenfelder, editor. Landscape Paradigms in Physics 
and Biology. Elsevier, Amsterdam, 1997. 

[3] D. J. Wales. Energy landscapes. Cambridge University 
Press, Cambridge, 2003. 

[4] G. Zarand, F. Pazmandi, K.F. Pal, and G.T. Zimanyi. 
Using hysteresis for optimization. Phys. Rev. Lett, 
89:150201, 2002. 

[5] K. F. Pal. Hysteretic optimization for the traveling sales- 
man problem. Physica A, 329:287-297, 2003. 

[6] K. F. Pal. Hysteretic optimization, faster and simpler. 
Physica A, 359:650-658, 2006. 

[7] K. F. Pal. Hysteretic optimization for the Sherrington- 
Kirkpatrick spin glass. Physica A, 367:261-268, 2006. 

[8] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi. Opti- 
mization by simulated annealing. Science, 220:671-680, 
1983. 

[9] S. Boettcher and A. G. Percus. Optimization with ex- 
tremal dynamics. Phys. Rev. Lett, 86:5211-5214, 2001. 

[10] D. Sherrington and S. Kirkpatrick. Solvable model of a 
spin-glass. Phys. Rev. Lett, 35:1792-1796, 1975. 

[11] S. F. Edwards and P. W. Anderson. Theory of spin 
glasses. J. Phys. F: Metal Phys., 5:965, 1975. 

[12] S. Zapperi, F. Colaiori, L. Dante, V. Basso, G. Durin, 
A. Magni, and M. J. Alava. Is demagnetization an effi- 
cient optimization method? J. Magn. Mat, E1009:272- 
276, 2004. 

[13] A.G. Percus, G. 1st rate, and C. Moore. Computational 
Complexity and Statistical Physics. Oxford University 
Press, New York, 2006. 

[14] G. Kotliar, P. W. Anderson, and D. L. Stein. One- 
dimensional spin-glass model with long-range random in- 
teractions. Phys. Rev. B, 27:R602, 1983. 

[15] D. S. Fisher and D. A. Huse. Ordered phase of short- 



range ising spin-glasses. Phys. Rev. Lett., 56:1601, 1985. 

[16] D. S. Fisher and D. A. Huse. Equilibrium behavior of the 
spin-glass ordered phase. Phys. Rev. B, 38:386, 1988. 

[17] H. G. Katzgraber and A. P. Young. Monte carlo studies 
of the one-dimensional ising spin glass with power-law 
interactions. Phys. Rev. B, 67:134410, 2003. 

[18] H. G. Katzgraber and A. P. Young. Probing the almeida- 
thouless line away from the mean-field model. Phys. Rev. 
B, 72:184416, 2005. 

[19] F. Pazmandi, G. Zarand, and G. T. Zimanyi. Self- 
organized criticality in the hysteresis of the Sherrington- 
Kirkpatrick model. Phys. Rev. Lett, 83:1034-1037, 1999. 

[20] H. H. Hoos and T. Stiitzle. Stochastic Local Search: 
Foundations and Applications. Morgan Kaufmann, San 
Francisco, 2004. 

[21] P. Salamon, P. Sibani, and R. Frost. Facts, Conjectures, 
and Improvements for Simulated Annealing. Society for 
Industrial & Applied Mathematics, 2002. 

[22] E. Marinari and G. Parisi. Simulated tempering: a new 
monte carlo scheme. Europhys. Lett., 19:451-458, 1992. 

[23] K. F. Pal. Hysteretic optimization, faster and simpler. 
Physica A, 360:525, 2006. 

[24] A. J. Bray, M. A. Moore, and A. P. Young. Lower criti- 
cal dimension of metallic vector spin-glasses. Phys. Rev. 
Lett, 56:2641-2644, 1986. 

[25] D. S. Fisher and D. A. Huse. Equilibrium behavior of 
the spin-glass ordered phase. Phys. Rev. B, 38:386-411, 
1988. 

[26] S. Boettcher and A. G. Percus. Nature's way of optimiz- 
ing. Artificial Intelligence, 119:275, 2000. 

[27] S. Boettcher. Extremal optimization for sherrington- 
kirkpatrick spin glasses. Euro. Phys. J. B, 46:501, 2005. 

[28] L. Viana and A. J. Bray. Phase diagrams for dilute spin- 
glasses. J. Phys. C: Solid State Phys., 18:3037, 1985. 



